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Abstract 

In  this  paper  two-dimensional  electromagnetic  scattering  problems  with  a  time-periodic 
incident  field  are  considered.  The  scatterer  is  a  perfect  conductor,  and  an  artificial  boundary 
condition  is  used.  The  large  time  behavior  of  solutions,  depending  on  (divergence-free) 
initial  conditions,  is  characterized.  It  turns  out  that  in  addition  to  the  expected  time- 
periodic  solution  the  limiting  solution  may  also  contain  a  spurious  stationary  field.  The 
source  of  the  stationary  field  is  explained  and  equations  describing  it  are  obtained.  Several 
avoidance  strategies  are  discussed,  and  numerical  comparisons  of  these  techiques  are  given. 
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1.  INTRODUCTION 


For  solving  scattering  problems  with  a  time-periodic  incident  field,  the  periodic  solution  is 
often  obtained  by  choosing  initial  conditions  and  time  marching  Maxwell’s  curl  equations 
to  a  periodic  state.  The  initial  conditions  must  satisfy  the  divergence  equations  but  are 
otherwise  arbitrary. 

When  this  time  marching  method  is  used  to  solve  scattering  problems  for  perfect  conductors 
it  is  found  that  the  computed  periodic  solutions  are  “contaminated”  by  spurious  stationary 
components.  For  example,  in  two  dimensional  transverse  magnetic  calculations  it  is  observed 
that  the  magnetic  field  is  affected  by  a  spurious  component;  in  three  dimensions  both  the 
magnetic  and  electric  fields  are  affected.  The  spurious  fields  can  be  relatively  large  and  must 
usually  be  removed  empirically  by  postprocessing,  for  example  by  peak-to-peak  averaging 
[2],  [3],  or  by  the  postprocessing  method  mentioned  in  section  4  below.  These  methods  are 
usually  effective  and  permit  quantities  such  as  radar  cross  sections  to  be  obtained  from  the 
resulting  numerical  data.  However,  postprocessing  introduces  additional  uncertainties  and 
it  is  worthwhile  to  try  to  obtain  accurate  numerical  solutions  directly. 

In  this  paper  we  wiU  explain  the  origin  of  the  spurious  fields  and  obtain  the  differential 
equations  and  boundary  conditions  which  characterize  them.  Based  on  these  results  we  will 
propose  some  techniques  which  enable  accurate  results  to  be  found  directly  and  efficiently 
without  postprocessing.  We  will  provide  numerical  comparisons  of  the  techniques. 

The  next  section  illustrates  the  occurrence  of  a  spurious  field  in  the  relatively  simple  setting 
of  a  two  dimensional  transverse  magnetic  computation.  Then  we  explain  the  source  of  this 
field  and  derive  the  equations  which  define  it.  Following  that  we  discuss  some  avoidance 
strategies  and  illustrate  their  properties.  The  remainder  of  the  paper  extends  the  results  to 
other  problems  and  briefly  mentions  some  mathematical  issues  which  are  analysed  in  detail 
in  [4]. 


2.  THE  MAGNETIC  OFFSET 


In  a  two  dimensional  transverse  magnetic  problem  it  is  the  magnetic  field  which  is  affected 
by  the  spurious  stationary  component.  To  illustrate  the  spurious  field  we  consider  the 
following  transverse  magnetic  scattering  problem:  in  fi. 


with  boundary  conditions 


dE 

=  curlH 

at 

(1) 

(2) 

div  H  =  0 

(3) 

— E,-  on  Fs 

(4) 

cyuH  X  n  =  0  on  To 

(5) 

E 

E 
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Figure  1:  The  time  history  of  (solid  line)  and  Hy  (dashed  line). 


and  initial  conditions 

E(x,  0)  =  Eo(x),  H(x,  0)  =  Ho(x).  (6) 

In  these  equations,  E  and  H  denote  the  scattered  field  variables,  E^-  is  the  incident  field, 
Ts  denotes  the  boundary  of  the  (perfect  conductor)  scatterer  and  To  is  the  artificial  outer 
boundary,  fl  denotes  the  domain  between  the  inner  and  outer  boundaries,  n  is  the  outer 
normal  to  and  x  =  (x,  y).  We  assume  that  E^-  can  be  written  in  the  form 

Ei(x,«)  =  e-*Ei(x), 

and  we  look  for  the  time-periodic  solutions  of  (l)-(6)  in  a  similar  form.  For  simplicity,  we 
are  assuming  a  first  order  radiation  condition  on  the  outer  boundary  .  The  results  are  not 
materially  affected  if  higher  order  radiation  conditions  are  used. 

To  illustrate  the  magnetic  offset  we  take  a  square  scatterer  with  sides  aligned  with  the 
coordinate  axes  and  electric  size  1  based  on  half  the  side  length  of  the  square.  The  scatterer 
is  illuminated  from  the  (1,0)  direction  by  a  field  of  the  form  The  outer 

boundary  is  ^  wavelengths  (4.5  body  lengths)  from  the  center  of  the  scatterer.  This 
problem  was  marched  to  a  periodic  state,  starting  from  zero  initial  conditions,  using  the 
standard  uniform  mesh  “FDTD”  scheme  [1]  with  20x  points  per  wavelength  and  with  a 
timestep  of  l/y/2  the  maximum  size  allowed  by  the  stability  condition. 

In  Figure  1  we  show  a  time  history  of  Hj,  and  Hj,  sampled  at  the  points  (-:^,|57)  and 
(— Mr  ’  x)’  respectively  (measured  from  the  center  of  the  body  in  wavelengths).  While  the 
records  appear  to  approach  periodic  functions,  they  are  not  symmetric  about  the  horizontal 
axis  and  consequently  they  cannot  be  of  the  desired  form  e’"*H(x).  A  similar  situation 
holds  at  the  other  mesh  points,  the  amount  of  the  offset  being  a  function  of  position.  This 
function  is  spurious  and  is  not  part  of  the  solution  to  the  periodically  forced  scattering 
problem.  We  wish  to  characterize  it  with  a  view  to  eliminating  it  altogether. 

In  contrast  to  Figure  1,  Figure  2  shows  the  time  history  of  the  electric  field  at  the  point 
(—7,  7).  There  is  no  obvious  problem  with  this  solution  and  in  fact  it  is  close  to  the  exact 


2 


Figure  2:  The  time  history  of  E, 


solution  of  the  periodically  forced  difference  equations.  A  general  proof  of  this  is  given  in  [4]. 


3.  ORIGIN  OF  THE  MAGNETIC  OFFSET 

Let  e“^*E  and  denote  the  periodic  solutions  obtained  by  solving  the  periodically  forced 
differential  equations  without  the  initial  conditions.  We  shall  assume  that  the  solutions  E 
and  H  of  the  initial  value  problem  are  such  that  for  large  times, 

E(x,  t)  -  e*‘^*E(x),  H(x,  t)  ~  e'‘^H(x)  +  H*(x).  (7) 

These  assumptions  reflect  what  is  seen  computationally. 

Since  e®"*E  and  satisfy  the  same  equations  as  E  and  H  (excluding  the  initial  condi¬ 

tions)  we  obtain  by  subtraction  the  following  equations  for  H*: 


divH* 

curlH* 

an* 

dt 

H*  X  n|r. 


0 

0 

0 

0. 


(8) 

(9) 

(10) 

(11) 


To  obtain  a  boundary  condition  on  the  scatterer  we  integrate  Faraday’s  law  with  respect 
to  time,  use  the  initial  conditions  and  evaluate  the  normal  components  on  Fg: 


1 

[(x)  •  n  =  Ho(x)  •  n  -  -  curl  E(x,  r)  •  n  dr  on  Tg. 


Since  curlE  •  n  on  Fg  is  the  tangential  derivative  of  E  along  the  boundary,  it  is  completely 
determined  by  the  boundary  condition.  So 


1  /** 

H(x)-n  =  Ho(x)-n-l--/  curlEi(x,r) •  n< 

Jo 


/  1  -  \ 

=  Ho(x)  •  n  —  — curlEi(x)  *11)  +  ^ — curlE{(x)  •  n  on  Fg. 

V  tojpi  J  tUjX 

In  the  above  expression  it  is  only  the  second  term  which  has  the  required  periodic  behavior 
and  it  follows  that  the  first  term,  which  is  stationary,  must  correspond  to  the  normal 
boundary  condition  for  H*.  Thus  we  are  led  to 

H*  •  n|r,  =  (Ho  •  n  -  ^curlE,-  •  n)|r,.  (12) 

To  interpret  this  boundary  condition,  assume  that  the  incident  magnetic  field  is  also  of  the 
form  e*‘^*H,(x)  such  that  Maxwell’s  equations  are  satisfied.  In  that  case  Faraday’s  law  gives 

e*‘^*H£(x)  =  T^e‘‘^*curlEi 

and  evaluating  this  at  t  =  0  and  using  it  to  replace  the  second  term  in  the  boundary  con¬ 
dition  it  follows  that  H*  •  n  is  the  normal  component  of  the  total  magnetic  field  at  t  =  0. 
Since  B  •  n  must  be  continuous  across  Fg  it  follows  that  there  is  a  magnetic  field  inside 
the  scatterer  if  the  right  side  of  (12)  is  nonzero.  In  addition  this  field  is  stationary  since 
E  is  zero  inside  a  perfect  conductor,  and  we  may  attribute  it  to  a  steady  current  flowing 
unhindered  in  the  conductor.  It  is  this  current  which  produces  the  magnetic  offset. 


4.  AVOIDING  THE  MAGNETIC  OFFSET 

In  this  section  we  wiU  discuss  four  different  possibilities  for  obtaining  the  correct  magnetic 
field  e®‘^*H(x).  The  first  approach  is  a  simple  postprocessing  technique  which  does  not 
use  the  results  of  the  previous  sections  except  for  the  assumptions  (7).  The  other  three 
approaches  do  rely  on  the  equations  obtained  in  the  previous  sections. 

Method  1:  Postprocessing.  The  postprocessing  approach  consists  simply  of  fitting  a 
function  of  the  form  A  -I-  to  the  computed  H  at  each  mesh  point,  and  then  discarding 
the  stationary  part  A.  A  decision  has  to  be  made  about  when  to  sample  H.  This  should  be 
at  reasonably  spaced  times  after  the  solution  is  believed  to  have  become  adequately  periodic 
(3-6  periods  is  often  sufficient). 

To  illustrate  the  stationary  part  we  took  a  square  scatterer  as  before  and  illuminated  it 
from  the  (1,0)  direction  with  a  field  Ei(x,t)  =  1000e*(‘^*“^‘’'-) .  The  problem  was  solved  (for 
the  imaginary  part  of  the  solution)  by  the  method  mentioned  above  in  section  2,  starting 
from  zero  initial  conditions.  We  calculated  (the  imaginary  part  of  )  the  stationary  part  by 
sampling  the  values  of  H  after  m,  m  —  ^  and  m  —  j  timesteps,  where  T  is  the  period  of 
the  incident  field  in  timesteps,  for  m  =  500  and  m  =  1000  (about  4  and  8  periods).  The 
values  of  the  stationary  parts  of  the  x-  and  y-components  of  H  on  lines  y  =  and  y  = 
respectively,  are  shown  in  figures  3  and  4.  Also  shown  are  the  theoretically  predicted  values 
of  the  stationary  field  H*  (the  approximate  solution  of  (8)-(12) ).  We  see  that  the  values  of 


Figure  3:  The  stationary  part  of  H^;.  Figure  4:  The  stationary  part  of  Hj,. 

Postprocessing  method  (zero  initial  data) 

500  timesteps  —  solid  line,  1000  timesteps  —  dashed  line,  H*  —  dotted  line. 


H*  are  quite  near  to  the  stationary  parts  of  H,  especially  for  longer  time.  The  stationary 
part  is  relatively  large  compared  to  the  amplitude  of  H,  which  is  near  one  for  most  of  ft. 

Method  2:  Solving  the  div-curl  system.  The  second  approach  to  elimination  of  the 
magnetic  offset  is  based  on  the  observation  from  (12)  that  if  the  initial  condition  for  H 

satisfies  Ho  •  n  =  ^ — curl  E,  •  n  on  Fg,  then  the  div-curl  system  for  H*  will  have  only  the 

tu/i 

zero  solution.  One  way  to  choose  a  suitable  Hq  is  as  a  solution  of  the  equations 


curl  Ho  =  0 

(13) 

div  Ho  =  0 

Ho  •  n  =  T^curl  E*-  •  n  on  Fg 
tuy, 

(14) 

(15) 

Ho  X  n  =  0  on  Fo. 

(16) 

To  illustrate  this  we  solved  Maxwell  equations  with  the  initial  conditions  Eq  =  0  and  Hq 
being  the  solution  of  (14)-(16).  Then  we  calculated  the  stationary  part  as  above.  The 
results  (the  values  of  the  stationary  parts  of  the  x-  and  t/-components  of  H  on  lines  y  = 
and  y  =  as  before)  are  shown  in  Figures  5  and  6. 

One  can  see  that  the  stationary  part  is  several  times  smaller  than  in  the  postprocessing 
method;  in  fact  it  is  probably  caused  by  the  discretization  error  in  time.  However,  even 
though  there  is  a  natural  way  to  formulate  the  div-curl  system  in  the  FDTD  framework, 
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Figure  5:  The  stationary  part  of  H,;.  Figure  6:  The  stationary  part  of  Hy. 

Ho  is  a  solution  of  the  div-curl  problem  (14)-(16). 

500  timesteps  —  solid  line,  1000  timesteps  —  dashed  line 


solving  it  can  be  quite  expensive.  Also,  up  to  round-off  error  the  solution  is  similar  to  that 
of  the  postprocessing  method  after  discarding  the  stationary  part.  So  it  seems  that  this 
method  is  not  worth  while. 


Method  3:  Extension  of  Ej.  Another  way  to  satisfy  the  condition  Ho-n  =  t— curlEj-n 

lOJH 

on  Ts  is  to  use  the  initial  conditions 


Eo(x)  =  -<^(x)E,(x),  Ho(x)  =  ^curl(^(x)Ei(x)),  (17) 

where  is  a  smooth  real- valued  function  such  that  <^  =  1  in  a  neighborhood  of  Fs,  and  ^  =  0 
in  a  neighborhood  of  Fq.  A  possible  advantage  of  this  scheme  is  that  since  the  initial  values 
and  the  boundary  data  match  at  t  =  0,  the  solution  should  be  smoother  in  time  for  small 
times.  We  would  then  expect  that  discretizations  should  have  a  better  rate  of  convergence 
as  the  mesh  size  approaches  zero. 

We  used  for  0  a  product  of  continuously  differentiable  cubic  splines  in  x  and  in  y.  The 
resulting  “stationary  field”  is  shown  on  Figures  7  and  8.  We  see  that  it  is  much  smaller  than 
in  the  previous  case,  especially  for  longer  time  (40-80  times  smaller  after  1000  timesteps). 

Method  4:  Modified  incident  field.  This  method  consists  of  giving  zero  initial  data 
and  making  an  adjustment  to  the  incident  field.  It  is  possible  to  change  this  field  for  a 
small  interval  of  time  so  that  the  spurious  field  will  be  zero.  One  possibility  is  to  define  the 
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Figure  7:  The  stationary  part  of  Figure  8:  The  stationary  part  of  H^. 

The  smooth  initial  data  (17). 

500  timesteps  —  solid  line,  1000  timesteps  —  dashed  line 


incident  field  to  be 


-€-^E,(x),  0<t<— , 
E,(x,i)  =  {  „ 


Then  at  time  to  =  on  Fg  we  have 
.  1  /•*»  t 


H(x,to)  •  n  =  -  /  curl  E,(x)  •  n  =  7-^e*“*“curl  E,(x)  •  n,  (19) 

Jo  to  tUfi  ^  ' 

so  that  H(x,to)  satisfies  the  correct  condition  at  time  to  =  Again  the  initial  and  the 
boundary  conditions  match,  so  the  solution  is  smoother.  Since  the  incident  field  reaches  the 
correct  value  after  the  first  period,  the  additional  work  required  is  roughly  equal  to  solving 
the  equations  during  one  period. 

The  results  are  shown  in  Figures  9  and  10. 

For  an  incident  field  in  the  form  it  might  be  preferable  to  adjust  the  field  as 

follows: 

f  0,  wt  <  k  •  X, 


Ei(x,  t)  - 


0  <  -  k  •  X  <  27r, 

ZTT 


27r  <  —  k  *  X. 
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Figure  9;  The  stationary  part  of  H^.  Figure  10:  The  stationary  part  of  Hj,. 

Continuous  adjustment  of  the  incident  field,  given  by  (18). 

500  timesteps  —  solid  line,  1000  timesteps  —  dashed  line 

Again  at  time  to  large  enough  so  that  the  incident  field  has  reached  the  correct  value  on  all 
Ts,  we  get  that  H(x,fo)  satisfies  the  condition  necessary  to  make  the  spurious  field  vanish. 
This  adjustment  makes  the  imaginary  part  of  E,-  continuously  differentiable,  so  we  should 
get  even  better  convergence  when  solving  for  the  imaginary  part  only.  The  results  are  shown 
in  Figures  11  and  12.  (We  used  m  =  600  instead  of  m  =  500,  because  after  4  periods  the 
adjusted  incident  field  has  barely  reached  its  correct  value  near  the  outer  boundary.)  Note 
that  the  stationary  parts  are  the  smallest  compared  with  the  other  methods. 

To  compare  the  smoothness  of  the  solutions  obtained  by  the  different  methods  we  present 
here  the  time  history  of  at  the  point  (-i,  ^).  This  point,  as  any  point  near  the  x-axis, 
is  special  in  the  sense  that  the  amplitude  of  Hj;  is  very  small,  so  effects  of  nonsmoothness 
stand  out  more  clearly.  Figure  13  shows  the  solutions  with  zero  initial  data  (method  1) 
and  in  the  case  of  “correct”  smooth  initial  data  (17)  (method  3).  The  graph  of  the  solution 
obtained  by  first  solving  the  div-curl  system  (method  2)  looks  the  same  as  of  the  solution 
obtained  by  method  1,  only  shifted  by  a  constant  (the  size  of  the  shift  is  determined  by 
the  value  of  the  solution  of  the  div-curl  problem  at  given  point),  so  we  chose  not  to  present 
it  here.  The  solutions  exhibit  very  large  oscillations  in  the  beginning,  They  correspond 
to  a  large  error  in  approximating  the  discontinuous  wavefront.  As  the  front  moves  away, 
the  solution  becames  smoother.  Then  the  oscillations  reappear  approximately  at  the  time 
when  the  (reflected)  wavefront  arrives  back  from  the  corners  of  the  artificial  boundary  (the 
reflections  from  the  smooth  part  of  the  boundary  are  probably  too  small  to  be  seen).  Figure 


Figure  11:  The  stationary  part  of  H^;.  Figure  12:  The  stationary  part  of  Hj,. 

Smooth  adjustment  of  the  incident  field,  given  by  (20). 

600  timesteps  —  solid  line,  1000  timesteps  —  dashed  line 


14  shows  the  solutions  with  zero  initial  data  and  altered  boundary  conditions  (method  4, 
boundary  conditions  (18)  and  (20)).  We  see  that  the  solution  looks  much  smoother  than 
the  solutions  obtained  with  the  other  methods,  so  one  might  expect  to  have  smaller  dis¬ 
cretization  error  than  in  the  other  cases. 
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Figure  13:  The  time  history  of  H^;  at  the  point  in  the  case  of  zero  initial  data 

(solid  line)  and  the  smooth  initial  data  (17)  (dashed  line). 


Figure  14:  The  time  history  of  Hj;  at  the  point  (—7,  4^)  in  the  case  of  continuous  ad¬ 
justment  of  the  boundary  conditions,  given  by  (18)  (solid  line),  and  in  the  case  of  smooth 
adjustment  (20)  (dashed  line),  both  with  zero  initial  data. 


5.  OTHER  PROBLEMS 

Consider  the  transverse  electric  problem:  in  0, 

'‘It  = 
div  E  =  0 

with  boundary  conditions 

E  X  n  =  — Ej-  X  n  on  Fg 

E  X  n  -f  c/xH  =  0  on  Fq 
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and  initial  conditions 


E(x,  0)  =  Eo(x),  H(x,  0)  =  Ho(x). 

In  this  case  there  is  an  electric  offset.  The  magnetic  offset  is  zero  because  H  •  n  =  0  on  Fg 
holds  automatically. 

Assuming  that  for  large  times 

E(x,  t)  ~  e*‘"*E(x)  +  E*(x),  H(x,  t)  ~  e~‘H(x), 
we  obtain  as  before  by  subtraction  the  following  equations  for  E*: 

dE* 

dt  ~  ^ 
curl  E*  =  0 

divE*  =  0 

E*  X  n  =  0  on  Fo  and  on  Fs 

A  system  such  as  this  is  determined  up  to  a  one  parameter  family  of  solutions.  In  this  case 
the  parameter  is  the  total  surface  charge  on  the  body.  This  quantity  should  be  zero  to  avoid 
a  spurious  addition  to  the  electric  field.  It  is  computed  as  follows: 

f  E-n=fEo -11  +  - f  [  curl  H  •  n  =  /  Eq  •  n. 

Jn  e  Jo  Jr,  Jr, 

It  follows  that  the  total  charge  is  constant  in  time  and  it  is  this  charge,  determined  by 
the  initial  condition,  which  is  responsible  for  the  electric  offset  field.  The  simplest  way  to 
eliminate  the  electric  offset  is  to  put  Eq  =  0.  The  system  above,  supplemented  with  this 
condition  of  zero  surface  charge  then  has  only  the  zero  solution. 

For  the  three  dimensional  case  the  calculations  above  remain  valid  so  both  the  electric  and 
magnetic  fields  may  have  spurious  components.  Our  equations  for  E*  and  H*  continue  to 
hold  and  so  should  the  remedies  although  we  have  not  checked  that  explicitly. 

It  has  not  been  possible  to  completely  rule  out  the  possibility  of  additional  (possibly  time 
dependent)  spurious  fields  in  three  dimensions  although  we  do  not  expect  that  such  fields 
exist. 

6.  MATHEMATICAL  COMMENTS 

We  have  made  a  number  of  assumptions  in  order  to  obtain  the  equations  for  the  offset  fields. 
We  wish  to  mention  the  major  assumptions  here. 

First  there  is  the  question  of  whether  the  equations  (l)-(6)  admit  a  periodic  solution.  This 
is  equivalent  to  the  question  of  existence  for  Helmholtz  equation  in  a  bounded  domain  with 
the  appropriate  boundary  conditions.  We  have  assumed  here  that  periodic  solutions  exist. 


11 


A  second  assumption  is  that  the  solutions  of  the  initial  value  problems  with  periodic  bound¬ 
ary  conditions  approach  to  a  limiting  state  as  t  — *■  oo.  While  this  is  certainly  expected  for 
reasonably  shaped  scatterers,  it  should  be  proved  in  a  precise  treatment. 

The  third  assumption  is  the  one  we  justified  by  appeal  to  computational  experiments,  which 
is  that  in  the  transverse  magnetic  (electric)  case  the  electric  (magnetic)  field  has  no  offset 
so  that  E(x,t)  -  e®‘^*E(x)  — »■  0  as  t  oo  (H(x,t)  -  e*"*H(x)  0)  as  t  oo.  Since  in  the 

transverse  magnetic  (electric)  problem  E  (H)  is  one-dimensional,  this  is  really  a  question  of 
the  asymptotic  behavior  of  the  wave  equation  with  periodic  Dirichlet  (Neumann)  condition 
on  the  inner  boundary  and  the  absorbing  boundary  condition  on  the  outer  boundary.  This 
has  been  known  for  a  while.  For  three  dimensions  the  question  is  more  delicate,  because 
the  equations  are  coupled  by  divergence  conditions,  and  also  because  both  the  electric  and 
the  magnetic  field  may  have  spurious  components. 

In  our  paper  [4],  which  is  a  companion  to  this  one,  we  provide  a  rigorous  statement  of 
these  assumptions  and  the  appropriate  results  are  proved  correct  for  the  two  dimensional 
problems. 


7.  CONCLUSIONS 

We  have  shown  the  origin  of  the  spurious  magnetic  and  electric  offset  fields  which  are 
encountered  when  time  marching  Maxwell  equations  for  periodic  solutions.  Based  on  the 
equations  governing  these  fields  several  remedies  were  discussed.  We  found  that  the  best 
remedy  is  to  modify  the  incident  field,  gradually  allowing  it  to  reach  its  desired  form  in  the 
way  specified  in  section  4.  This  technique  reliably  reduces  the  magnetic  offset  to  a  negligible 
size  and  avoids  the  need  for  postprocessing  of  the  results  of  the  computation. 
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